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We use Monte Carlo simulations to study polymer melts consisting of fully flexible and moderately 
stiff chains in the bond fluctuation model at a volume fraction 0.5. In order to reduce the local 
density fluctuations, we test a pre-packing process for the preparation of the initial configurations of 
the polymer melts, before the excluded volume interaction is switched on completely. This process 
leads to a significantly faster decrease of the number of overlapping monomers on the lattice. This 
is useful for simulating very large systems, where the statistical properties of the model with a 
marginally incomplete elimination of excluded volume violations are the same as those of the model 
with strictly excluded volume. We find that the internal mean square end-to-end distance for 
moderately stiff chains in a melt can be very well described by a freely rotating chain model with a 
precise estimate of the bond-bond orientational correlation between two successive bond vectors in 
equilibrium. The plot of the probability distributions of the reduced end-to-end distance of chains 
of different stiffness also shows that the data collapse is excellent and described very well by the 
Gaussian distribution for ideal chains. However, while our results confirm the systematic deviations 
between Gaussian statistics for the chain structure factor S c {q) [minimum in the Kratky-plot] found 
by Wittmer et al. {EPL 77 56003 (2007).} for fully flexible chains in a melt, we show that for the 
available chain length these deviations are no longer visible, when the chain stiffness is included. 
The mean square bond length and the compressibility estimated from collective structure factors 
depend slightly on the stiffness of the chains. 

I. INTRODUCTION 

In the theoretical study of polymer physics 00 , com¬ 
puter simulations provide a powerful method to mimic 
the behavior of polymers covering the range from atomic 
to coarse-grained scales depending on the problems one 
is interested in Ml- The generic scaling properties of 
single linear and branched polymers in the bulk under 
various solvent conditions and polymer solutions at dif¬ 
ferent temperatures have been described quite well by 
simple coarse-grained lattice models, e.g. self-avoiding 
walks on the simple cubic lattice, and off-lattice models, 
e.g. bead-spring model, etc. A wide variety of computa¬ 
tional strategies has been employed to simulate and an¬ 
alyze these models, including conventional (Metropolis) 

Monte Carlo (MC) schemes with various types of moves, 
chain growth algorithms, parallel tempering, and molec¬ 
ular dynamics [3[. On the one hand, however, as the size 
and complexity of a system increases, detailed informa¬ 
tion at the atomic scale may be lost when employing low 
resolution coarse-graining representations. On the other 
hand, the cost of computing time may be too high if the 
system is described at high resolution. Therefore, more 
scientific effort has been devoted to developing an appro¬ 
priate coarse-grained model where the connection to the 
details on the atomic scale is not lost, but rather atom¬ 
istic details are suitably mapped on effective potentials 
on the coarser scale. Such models then can reproduce the 
global thermodynamic properties and the local mechan¬ 
ical and chemical properties such as the intermolecular 
forces between polymer chains Ml- Improving such 
approach further is still an active area of research. 

In this work we deal with linear polymer chains in a 
melt on the simple cubic lattice. Although coarse-grained 
lattice models neglect the chemical detail of a specific 


polymer chain and only keep chain connectivity (topol¬ 
ogy) and excluded volume, the universal behavior of poly¬ 
mers still remains the same in the thermodynamic limit 
(as the chain length N —> oo)[2.]. We consider for our 
simulations the bond fluctuation model (BFM) @,0- 
E3 where linear chains are located on the simple cubic 
lattice with bond constraints. The BFM has the advan¬ 
tages that the computational efficiency of lattice mod¬ 
els is kept and the behavior of polymers in a continuum 
space can be described approximately. The model thus 
introduces some local conformational flexibility while re¬ 
taining the computational efficiency of lattice models for 
implementing excluded volume interactions by enforcing 
a single occupation of each lattice vertex. Although much 
work using this model exists already, only the fully flex¬ 
ible limit of the model has been used exclusively. This 
limit does not suffice when one considers a possible map¬ 
ping of atomistic details to this model, which requires to 
include some description of chain stiffness as well. 

A review of recent BFM studies is given in Ref. Il8j . 
When the concentration of polymer solutions is above the 
overlap concentration denoted bye*, the excluded vol¬ 
ume interactions are screened 0 [lj| . The average inter¬ 
action between monomers finally should cancel in a poly¬ 
mer melt since every monomer is isotropically surrounded 
by other monomers belonging to the same chain or not 
according to Flory’s argument. Therefore the polymer 
chains in a melt behave as ideal chains, where the ex¬ 
cluded volume effect is no longer important. However, 
a careful investigation of an individual polymer chain in 
a polymer melt based on the BFM and the bead-spring 
model (BSM) shows that there are noticeable deviations 
from an ideal chain behavior 0, |2Cj. The deviations 
are due to the incompressibility constraint of the melt. 
Therefore, small scale-free corrections to the asymptotic 
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behavior for ideal chains exist and they are irrespective 
of the stiffness of the chains. In Refs. [l 8 j, [2(|, the au¬ 
thors have shown that the corrections must decay as the 
dimensionless (Ginzburg) parameter p* / p ~ 1 /(pfr^s 1 / 2 ) 
for 1 -C s <C N with the increase of the segments s 
and the effective bond length b e of the chain. Here p is 
the overall monomer density and p* ~ s/R(s) 3 the seg¬ 
mental overlap density related to the internal distance 
R(s) ~ b e s 1 / 2 . It is still unclear whether the predicted 
power-law deviations from unperturbed Gaussian behav¬ 
ior can be observed for the internal segment lengths of 
semiflexible chains in the range 1 -C g{£b) "C s -C N by 
varying the stiffness of chains, i.e. tuning the bending 
energy £&. 

A very important practical problem for the simulation 
of melts of very long polymer chains is the construction of 
a suitable initial configuration. A naive approach for cre¬ 
ating initial configurations of melts with long polymers, 
would be to switch off initially the excluded volume inter¬ 
actions between monomers (so that the polymers resem¬ 
ble Gaussian chains), arrange these chains randomly in 
a simulation box, and switch on gradually the excluded 
volume interactions at the final step. It has been demon¬ 
strated by simulations using the BSM in the continuum 
that this strategy is not feasible since it leads to chain 
deformations on short length scales [2l| • With increasing 
chain length this effect becomes more pronounced due to 
the fact that self-screening and correlation hole effects, 
inherent to flexible high polymer melts, are not properly 
accounted for by this procedure. This problem was over¬ 
come by reducing the density fluctuations of Gaussian 
chains via a “pre-packing” strategy, and then switching 
on the excluded volume interactions in a quasi-static way 
(slow push-off) HHHH- We test a similar strategy of sim¬ 
ulating polymer melts based on the lattice model BFM 
in this work to check whether it helps to generate nearly 
equilibrated initial configurations through such a proce¬ 
dure. 

The outline of the paper is as follows: Sec. ITT1 describes 
the model and the simulation technique discussing how 
to generate initial configurations of polymer chains in 
a melt to equilibrate the system. Sec. um presents the 
simulation results at different stages from the initial state 
to the equilibrium state of fully flexible and moderately 
stiff polymer chains in a melt. Finally our conclusions 
are summarized in Sec. hyi 


II. MODEL AND SIMULATION TECHNIQUES 

_In the standard bond fluctuation model (BFM) 0 0 - 

iia a flexible polymer chain with excluded volume in¬ 
teractions is described by a self-avoiding walk (SAW) 
chain of effective monomers on a simple cubic lattice 
(the lattice spacing is the unit of length). Each ef¬ 
fective monomer of such a SAW chain blocks all eight 
corners of an elementary cube of the simple cubic lat¬ 
tice from further occupation. Two successive monomers 


along a chain are connected by a bond vector b which is 
taken from the set {(± 2 , 0 , 0 ),(± 2 , ± 1 , 0 ), (± 2 ,± 1 ,± 1 ), 
(±2,±2,±1), (±3,0,0), (±3,±1,0)} including also all 
permutations. The bond length | b | is therefore in a 
range between 2 and VTO. There are in total 108 bond 
vectors and 87 different bond angles between two sequen¬ 
tial bonds along a chain serving as candidates for building 
the conformational structure of polymers. 

As the stiffness of the chains is considered, semiflexible 
chains of N monomers in a melt thus are described by 
SAW chains on the simple cubic lattice, with a bending 
potential 
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where £& is the bending energy (ej, = 0 for ordinary 
SAWs), and 6^+1 is the bond angle between the i th bond 
vector and the ( i ± l) th bond vector along a chain. k B T 
is of order unity throughout the whole paper. A sys¬ 
tematical study of single semiflexible chains based on the 
BFM [23| shows that due to bond vector fluctuations and 
lattice artifacts the initial decay of the bond-bond ori¬ 
entational correlation function deviates from the simple 
exponential decay for £& > 10. Therefore, one should be 
careful of using the BFM for studying rather stiff chains 
as has already been pointed out also in Ref. [24| . 

It is well understood that linear polymers in a melt 
can be described by SAW chains based on the BFM at 
a volume fraction (j) = 0.5 @, IH, E3- Since in the BFM, 
each effective monomer occupies one unit cell, containing 
8 lattice sites, the monomer density is therefore defined 
as p = 4>/8. A detailed description of how we prepare 
the initial configurations is given as follows: The initial 
configuration of polymer melts containing n c chains of N 
monomers in a box of size V = L 3 with periodic bound¬ 
ary conditions in all three directions are generated in the 
following way: 

• The linear dimension of a simple cubic lattice is 
set up as L = (8n c Nj <^) 1//3 with volume fraction 
cj) = 0.5. 

• At the 0 th step the first n c monomers of the chains 
are randomly put on the lattice sites without dou¬ 
ble overlapping. It can be easily done by divid¬ 
ing the box into n c blocks and putting a single 
monomer in each block at a randomly chosen posi¬ 
tion inside the block. 


• Polymer chains are built like non-reversal random 
walks (NRRWs) by adding one monomer at each 
step until all n c chains reaching the required chain 
length N — 1. (We define here the chain length as 
the number of bonds along the contour of a chain. 
The contour length of the chain then is L = (N — 
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1 )£b where lb = (| b 2 I) 1 / 2 is the root-mean-square 
bond length). Based on the BFM there are 108 
possibilities to place the next monomer for each 
chain at the 1 st step, but at the following steps 
only one of 107 directions can be selected since an 
immediately reverse step is not allowed. At this 
stage, the excluded volume effect is switched off 
completely and n c NRRWs of (TV — l)-steps in the 
simulation box are generated. 

If the stiffness of the chains is considered, the prob¬ 
ability of placing the (i + l) th monomer connected 
to the i th monomer is proportional to exp(—£;,(1 — 
cos fli.i+l)). 

• Finally the excluded volume interactions between 
monomers are considered by applying the Local 
26 and slithering-snake moves to relax the poly¬ 
mer chains and push off those monomers blocking 
the same lattice site until all chains satisfy self- 
and mutual-avoidance. This “steepest descent” ap¬ 
proach only works for lattice chains since excluded 
volume constraints can be checked very efficiently, 
but would not work for continuum chains. 


Due to the large local density fluctuations of randomly 
overlapping NRRW chains in a melt, a pre-packing pro¬ 
cedure to reduce the chain deformations is also tested 
before the excluded volume interactions are switched on. 
This pre-packing procedure has been successfully per¬ 
formed in the simulations of equilibrating very long poly¬ 
mer chains in melts up to n c = 1000 and N = 1000 based 
on the bead-spring model in the continuum [22? . In this 
procedure all NRRW chains still keep their own struc¬ 
tures as rigid bodies, but are rearranged by Monte Carlo 
(MC) trial moves. The cost function which describes the 
average fluctuations of the particle number n(er) within a 
sphere of radius a with its center located at any monomer 
i is defined by 


U(cr) = < n(a) 2 > — < n(a) > 2 


( 2 ) 
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( 3 ) 

where N tot = Nn c is the total number of monomers in 
a melt, rii{a) is the number of monomers inside the i th 
sphere of radius cr, r,;j =| r t — fj | is the distance between 
monomers i and j irrespective of the chain connectiv¬ 
ity, and H{x) is the Heaviside step function. A trial 
move is accepted if the cost function becomes smaller, 
namely, the local density fluctuation is minimized in this 
process. Two types of MC moves, translation and pivot¬ 
like moves, applied to the center of mass of an individ¬ 
ual chain are considered as follows: (a) In the transla¬ 
tion move, an individual chain (or its center of mass) 




0 , for x < 0 
x > 0 


is randomly translated by a vector chosen from the set 
{(±1, 0,0), (±1, ±1, 0), (±1, ±1, ±1)} including all per¬ 
mutations. (b) In the pivot-like move, a whole individual 
chain is transformed by randomly choosing one of the 47 
symmetric operators including rotations by 90° or 180° 
around a random axis and reflections through its center 
of mass, and inversion at its center of mass. The cen¬ 
ter of mass serves as a pivot point here. Note that all 
monomers are only allowed to sit on the lattice sites, so 
a chain would need to be shifted slightly after the pivot¬ 
like move. However, the deviation of its center of mass 
is less than 1 lattice spacing. 

For equilibrating the system and measuring interesting 
physical quantities in the equilibrium runs three types of 
MC moves, local 26 moves [25j], slithering-snake moves, 
and pivot moves, are used in our simulations and briefly 
described as follows, 

• Local 26 (“L26”) move: a monomer is chosen 
randomly to move to one of the 26 nearest and 
next nearest neighbor sites surrounding it, i.e, ran¬ 
domly translated by a vector chosen from the set 
{(±1, 0,0), (±1, ±1, 0), (±1, ±1, ±1)} including all 
permutations. The bond crossing is allowed during 
the move. This is different from the traditional lo¬ 
cal 6 (“L6”) move where a monomer of the chain is 
chosen at random and moved to the nearest neigh¬ 
bor sites in the six lattice directions randomly. 

• Slithering-snake move: an end monomer is removed 
randomly and connected to the other end of the 
same linear chain by a bond vector randomly cho¬ 
sen from the set of 108 allowed bond vectors. 

• Pivot move: a monomer is chosen randomly from a 
linear chain as a pivot point, and the short part of 
the linear chain is transformed by randomly select¬ 
ing one of the 48 symmetry operations (no change; 
rotations by 90° and 180°; reflections and inver¬ 
sions). 

Of course, any attempted Monte Carlo move is accepted 
only if it does not violate constraints of our physical 
systems, such as excluded volume and bond length con¬ 
straints. 


III. SIMULATION RESULTS 

We first wish to test whether applying the pre-packing 
process to rearrange the NRRW chains based on the BFM 
before the excluded volume interaction is switched on 
helps to prepare the nearly equilibrated initial config¬ 
urations or not. We restrict the size of the lattice to 
be L 3 = 128 3 in our simulations. The total number 
of effective monomers based on the BFM is therefore 
-Ntot = n c N = 131072. Varying n c and N but keep¬ 
ing Ntot fixed the conformational properties of polymer 
melts including fully flexible chains (£& = 0.0) and mod¬ 
erately stiff chains [eb = 1-0, 2.0, and 3.0) are studied. 
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FIG. 1: Time series of local density fluctuation ((n(cr) 2 ) — 
(n(<r)) 2 )/(n(cr)) for a = 12. Data are for n c = 256 NRRW 
chains of IV = 512 monomers in a box of size L 3 = 128 3 . 
Several bending energies Eb are chosen, as indicated. 


Results of simulating polymer melts from preparing the 
initial configuration to the analysis of equilibrated chain 
structures in a melt are discussed in the next Subsection 
III A. 


A. Local density fluctuations 

In Fig. [U we show the typical time series of the lo¬ 
cal density fluctuation ((n(cr) 2 } — (n(a)) 2 ) / (n(a)) with 
cr = 12 for n c = 256 NRRW chains of N = 512 monomers 
in a box of size L 3 = 128 3 , and for the bending energies 
Eb = 0 (fully flexible), 1.0, 2.0, and 3.0. In the pre¬ 
packing process, the radius a is first set to 12 « 4 £b 
in the unit of lattice spacings until all curves describ¬ 
ing the local density fluctuation reach a plateau, and 
then reduced to 6 ~ 2 ib until all curves reach another 
plateau again. The estimates of the collective structure 
factor S(q) for q = 2i t/L and the local density fluctua¬ 
tion ((n(a) 2 ) — (n(a)) 2 ) / (n(a)) for ex = 12 are of the same 
order of magnitude (see Sec. IIIIDD from the final config¬ 
urations generated by this two-step pre-packing process. 
Each Monte Carlo step contains n c translation moves and 
n c pivot-like moves such that the locations of all chains 
in the box are rearranged by the trial moves. Each curve 
shown in Fig. [T| represents the results averaged over 32 
independent realizations. We see that the local density 
fluctuation can be reduced by about an order of magni¬ 
tude after about 1000 Monte Carlo steps. 

It is interesting to know what we have gained by in¬ 
troducing the pre-packing procedure before equilibrating 
the system. Therefore, we check the process of push¬ 
ing away the monomers occupying the same unit cells 
(eight lattice sites are occupied by one effective monomer 
in BFM) by observing the time series of the number of 
overlapping monomers, N oveT , as the excluded volume 
effect is switched on. For simplicity, if a lattice site is 
occupied by s effective monomers, it contributes s — 1 


overlapping monomers to N ovel . The transition process 
from n c = 256 NRRW chains of N = 512 monomers 
to SAW chains is observed according to the decrease 
of IVover shown in Fig. [2] for eb = 0.0, 1.0, 2.0, and 
3.0. Two sets of data (8 samples for each) represent¬ 
ing the push-off process before (dashed curve in green 
color online) and after (solid curves in red color online) 
the pre-packing process is applied. Here one MC step is 
a sequence of N tot L26 moves, and Nn c slithering-snake 
moves. In the slithering-snake moves, an end-monomer of 
each chain is selected and attempted to move once until 
all chains have been chosen, and then the same proce¬ 
dure is repeated N times. A trial move is accepted if 
AAover = Aive e r w) - N.< 0 and exp(-A Ub)/k B T > r 
where r is a random number and r £ [0,1). The latter 
condition is the Metropolis criterion preserving the stiff¬ 
ness of the chains. If only L26 moves are applied in this 
push-off process, at the beginning it is very efficient to 
push off overlapping monomers, but then some chains are 
trapped in a state containing knots, or a state blocked 
by other chains. The slithering-snake moves therefore in¬ 
deed play an important role at the intermediate stages 
for the chains escaping from a trapped state. From our 
observations in Fig. [2] it seems that applying an addi¬ 
tional pre-packing process does not speed up the push-off 
process. However, it is still necessary to investigate the 
conformation of chains in a melt carefully. 

In addition, it is interesting to note that pre-packing 
leads to a significantly faster decrease of A over with time 
during intermediate stages of the process. This feature 
may be of interest in cases where one does not require 
that strictly N ovel = 0 in the initial stage of the av¬ 
eraging. Consider e.g. a variant of the model where 
monomers overlap is not strictly forbidden (which corre¬ 
sponds to an infinitely high repulsive energy) but only 
leads to a large but finite energy penalty, e.g. E ovel = 
20 ksT. The probability to accept a move which increases 
the overlap then is of order 10 -9 , i.e. essentially negligi¬ 
ble, and for long chains the statistical properties of the 
model are essentially the same as for the model with strict 
excluded volume. For such a model the overlap energy 
per monomer reduces after about 10 4 Monte Carlo steps 
from an initial value (of the order of 20 ksT) by three or¬ 
ders of magnitude, i.e. 0.02 ksT, and this small number 
is reached 3 to 4 times faster with pre-packing than with¬ 
out it. For very large systems (and complicated models, 
leading to slow simulation algorithms) it may hence be 
an acceptable compromise to allow a small fraction (10~ 3 
or less) of overlap in the system. 

In Table Q] we compare the average values of local den¬ 
sity fluctuation ((n(a) 2 ) — (n(a)) 2 )/(n(a)) {Eq. J3])} for 

a = 12 (since g m i n = 2tt/L ss Kut/ 3 ) over 32 con¬ 
figures at different stages from the preparation of the 
initial configurations to the final configurations at the 
equilibrium state. Two ways of preparing the initial 
configurations (Sec. |TT| are compared as follows: (1) 
Randomly generating weakly-interactive NRRWs (NR- 
RWs), applying the pre-packing process (Pre-packing, 
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(a) (b) 






FIG. 2: Time series of the total number of overlapping monomers N ovel for polymer melts consisting of n c = 256 chains of 
N = 512 monomers, and for e b = 0.0 (a), 1.0 (b), 2.0 (c), and 3.0 (d). Two sets of data for the simulations with (solid curves 
in red color online) and without (dashed curves in green color online) the pre-packing process, as indicated. 


NRRWs), and switching on the excluded volume effect 
(Pre-packing+EV, SAWs). (2) Randomly generating 
weakly-interactive NRRWs (NRRWs) and switching on 
the excluded volume effect (EV, SAWs). We see that 
for weakly-interactive NRRWs the local density fluctua¬ 
tion decreases as the bending energy e b increases. Since 
the bond angles tend to become smaller, the distance be¬ 
tween non-bonded monomers apart from each other along 
the chain becomes longer. After the pre-packing process 
the fluctuation reduces by about 95%. Once the excluded 
volume interactions between monomers are switched on, 
the fluctuation remains almost the same for all cases. 


distance, (i? 2 (s)}, 

l N ~ s ^ ^ ] \ 

N _ s E ~ r ".i+s) 2 ) , 

3 = 1 / 

(4) 

where s is the chemical distance between the j th 
monomer and the (j + s) th monomer along the identical 
chain. The theoretical prediction of the internal mean 
square end-to-end distance for polymer melts consisting 
of semiflexible chains in the absence of excluded volume 
effect described by a freely rotating chain model is [l] 


(R 2 (s)) = 


■E 

n —1 


B. Mean square internal end-to-end distance and 
bond-bond orientational correlation function 


(R 2 (s))=s£ 2 

1 + (cos 6) 

2 (cos0)(l — ((cos0)) s 


1 — (cos 9) 

s(l — (cos0)) 2 

5 

(5) 


with 




f b = (?) 

and 

(cos 9) = (bi ■ bi+i)/£ 2 

(6) 


For understanding the connectivity and correlation be¬ 
tween monomers the conformations of linear chains of 
contour length L = (TV — l)£ b in a melt are normally de¬ 
scribed by the average mean square internal end-to-end 


where t b is the root-mean-square bond length. In the 
limit TV —> 00 , the bond-bond orientational correla¬ 
tion function in the absence of excluded volume effects 
therefore decays exponentially as a function of chemi¬ 
cal distance s between any two bonds along a linear 
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TABLE I: Values of local density fluctuation (( n 2 ) — (n) 2 )/(n) for a = 12. Data are for N = 512, and for £& = 0.0, 1.0, 2.0, 
3.0, averaged over 32 configurations at different stages. 


£b 

0.0 

1.0 

2.0 

3.0 

NRRWs 

51.06 

30.96 

20.89 

14.57 

Pre-Packing, NRRWs 

2.92 

1.56 

1.04 

0.86 

Pre-Packing+EV, SAWs 

0.28 

0.28 

0.28 

0.28 

EV, SAWs 

0.27 

0.27 

0.28 

0.27 

Equilibrium, SAWs 

0.28 

0.28 

0.27 

0.27 


chain 0,0, 

{bi ■ b i+s ) = £l(cos9(s)) = £l(cos9) s = i\ exp(-s4/4), 

(7) 

where £ p is the so-called persistence length which can be 
extracted from the initial decay of (cos 9(s)}. Equiva¬ 
lently, one can calculate the persistence length from 

£ P ,e/h = -1/ ln((cos 9)) (8) 

here instead of £ p we use £ p ,0 to distinguish between these 
two measurements. Replacing s by TV — 1 in Eq. ([5]) it 
gives the asymptotic behavior of the mean square end- 
to-end distance of a FRC equivalent to the behavior of a 
freely jointed chain 

(Rj) = C^N-iyl with Coo = [^21] (9) 

= n K i% : = 2 i p L (10) 

where Coo is called Flory’s characteristic ratio 

2 l p is the Kuhn length, and nx is the number of Kuhn 

segments. 

We compare the estimates of ( R 2 (s)) obtained from the 
initial configurations of SAW chains generated by push¬ 
ing off monomers occupying the same unit cell before and 
after the pre-packing process to that from those configu¬ 
rations in equilibrium as shown in Fig. [3] The two ways 
of preparing the initial configurations with and without 
the pre-packing process are denoted by “EV” and “Pre- 
packing-1-EV” , respectively in the figure. Results are ob¬ 
tained by taking the average over 32 initial configurations 
of SAW chains for each case. We see that there is only 
a slight discrepancy in the estimates of ( R 2 (s )) obtained 
from the initial configurations and from the equilibrated 
configurations. The discrepancy becomes more promi¬ 
nent for £b = 3.0. We see that there is no significant 
advantage of preparing the initial configurations of lat¬ 
tice chains in a melt through the pre-packing process as 
that was observed for the continuum chains. 

In the equilibration process, each Monte Carlo step 
contains n c N L26 moves where each monomer is selected 
once to move, n c slithering-snake moves where each chain 
is selected once to move, and n c pivot moves where each 
chain is also selected once to move. There are about 
10 6 independent configurations for each measurement in 



FIG. 3: Rescaled internal mean square end-to-end distance, 
{R 2 (s)) / (si 2 ), plotted as a function of s for polymer melts 
containing n c = 256 chains of N = 512 monomers, and for 
Eb = 0.0, 1.0, 2.0, and 3.0. Data are for the estimates ob¬ 
tained from the initial configurations of SAWs generated by 
the two methods, ”Pre-packing+EV”, and ”EV”, and from 
the equilibrated configurations. 


equilibrium. Figure [I] represents the well-equilibrated 
data for polymer melts in a box of size L 3 = 128 3 . Re¬ 
sults for fully flexible chains of chain sizes N = 32, 64, 
128, and 512 shown in Fig. QJi are in perfect ag reement 
with the theoretical prediction given in Ref. [18l . |25| , 

R 2 (s) = &e s (l — Cs/yfs) for 1 « s « N (11) 

where b e = lim s _,. 00 (l? 2 (s)/s) 1 / 2 is the “effective bond 
length” (which is different from the root-mean-square 
bond length 4 between monomers along a chain), and 
c s = y/24j 7r 3 /(p6 3 ). Results for polymer melts of N = 
128, 256, and 512 monomers are shown in Fig. 0 Jd for 
£b = 0.0, 1.0, 2.0, and 3.0. The theoretical prediction for 
a freely rotating chain (FRC), Eq. dSJ), is also included for 
comparison. For £& = 0.0 we see a strong deviation from 
the theoretical prediction since chains tend to swell due to 
the excluded volume effect that causes the 1 / y/s correc¬ 
tion in Eq. (fill) . This is different from the results [U 
obtained using the bead-spring model in the continuum, 
where on short and intermediate length scales Eq. © 
for a FRC overestimates the internal distances for fully 
flexible chains, while on a large length scale it gives the 
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FIG. 4: (a) (R 2 (s))/s plotted versus s for polymer melts consisting of n c N = 131072 monomers, and for Eb = 0.0. (b) 
(R 2 (s))/ (s£ 2 ) plotted versus s for polymer melts of IV = 128, 256 and 512 monomers, and for various values of bending energy 
£(,, as indicated. The theoretical predictions, Eq GD with b e = 3.244 and c 3 = 0.412 in (a), and Eq. © in (b), are shown by 
curves for comparison. 


(a) (b) 




FIG. 5: (a) Semi-log plot of the bond-bond orientational correlation function (cos d(s)) vs. s. (b) Same data as shown in (a), 
but on a log-log scale. Data are for polymer melts consisting of semiflexible chains of size N = 512, and for Eb = 0.0, 1.0, 2.0 
and 3.0, as indicated. In (a) the straight lines indicate the fits of the initial decay, (cos 9(s)) oc exp (—s£b/£p) (Eq. (7)), for 
Eb > 0. For Eb ■— 0 the straight line is described by exp(— slb/£ P ,e)- Values of i v /tb and l Pt e/£b are listed in Table III. In (b) 
the straight line indicates the power law, (cos d(s)) oc s~ 3 ^ 2 . 


accurate prediction. As the bending energy £{, increases, 
we see in Fig. 0 )d that the deviation from the prediction 
for a FRC reduces. 

We also estimate the root-mean-square bond length 
£b = (6 2 ) 1 / 2 , end-to-end distance (i? 2 ), and radius of gy¬ 
ration (i? 2 ). Results for various values of £b and TV are 
listed in Tabic HT1 and IIIII We see that (i? 2 )/(i? 2 ) ~ 6 
which is predicted for ideal chains. Our results of (i? 2 ), 
(i? 2 ), and £b for £b = 0.0 and for various values of TV are 
in perfect agreement with the results given in Ref. J2j| al¬ 
though the total number of monomers of polymer chains 
in a melt and the lattice size which we have chosen for our 
simulations both are eight times smaller. The persistence 
lengths for polymer melts consisting of semiflexible chains 
determined by fitting the exponential decay, Eq. 0, and 
the correlation between two neighboring bonds, Eq. 0, 


are shown in Fig. [5ji and also listed in Table EH includ¬ 
ing Flory’s characteristic ratio Coo using Eq. @. Note 
that the asymptotic decay of (cos 9(s)) with s is not ex¬ 
ponential as predicted by Eq. 0, but rather a power law 
decay, (cos 0(s)) oc s -3 / 2 for s* < s -C TV, is expected, 
due to excluded volume effects [25|, HI| . 

Therefore only the initial decay of (cos d(s)) with s 
(s < s*) is meaningful for the estimation of the per¬ 
sistence length. In Fig. [5]o we plot (cos0(s)} versus s 
on a log-log scale. The power law s -3 / 2 is also shown 
for comparison. We see that the range over which the 
power law behavior holds (s > s*) shrinks with increas¬ 
ing chain stiffness, and the data deviate from the straight 
line showing that the finite-size effect sets in. Results 
for various values of £b and for chain sizes TV = 128 
and N — 512 are listed in Table ED The relationship 
tp/lb ~ £ P ,e/h ~ Coo/2 holds for semiflexible chains 













TABLE II: Estimates of root-mean-square end-to-end distance {R 2 ) 1 ^ 2 , gyration radius ( R 2 ) 1 ^ 2 , and bond length lb = (6 2 ) 1 / 2 
for polymer melts containing n c fully flexible chains of N monomers (sb = 0). 


N 

n c 

{Rl) 1/2 

V6{R 2 ) 1/2 

lb 

32 4096 

17.07 

17.13 

2.63 

64 2048 

24.74 

24.75 

2.63 

128 

1024 

35.54 

35.51 

2.63 

256 

512 

50.76 

50.73 

2.63 

512 

256 

72.17 

72.20 

2.63 


(£b > 0) as predicted in Eqs. © and (fTUl) . The root- 
mean-square bond length lb ~ 2.63 for polymer melts, 
shown in Table hd is smaller than lb ~ 2.72 for linear 
flexible chains in dilute solution based on the BFM. The 
reason is that the crowded environment of polymer chains 
prevents the bond vector between two effective monomers 
going through several unit cells, i.e. occupying large vol¬ 
ume. We also see that lb depends only weakly on the 
stiffness of the chains as shown in Table ED 


C. Probability distributions of R e and R g 

The conformational behavior of polymer chains of size 
IV in a melt can also be described by the probability dis¬ 
tributions of end-to-end distance R e and gyration radius 
Rg, Pn{.R e) and P/v(P g ), respectively. The probability 
distribution of R e for ideal chains is simply a Gaussian 
distribution, 


p,A=f = pf-^%). m 

V e) \2nNl 2 J 2 Nl 2 J 

There exists exact theoretical prediction of the prob¬ 
ability distribution of R g for ideal chains [ 29 M 3 TI ] . but 
they are complicated to evaluate. However, it has been 
checked fill l23l l32j that the same formula suggested by 
Lhuillier [331 ] for polymer chains under good solvent con¬ 
ditions in d-dimensions is still a good approximation for 
ideal chains, i.e., 


Pn{.R 9 ) ~ exp 



( 13 ) 


where a± and 02 are (non-universal) constants, and the 
exponents a and <5 are linked to the space dimension d 
and the Flory exponent v by a = (vd — l)^ 1 and 5 = 
(1 — v)- 1 . Here (l + cc) is the des Cloizeaux exponent [H] 
for the osmotic pressure of a semidilute polymer solution, 
and S is the Fisher exponent [35| characterizing the end- 
to-end distance distribution. 

Numerically, the probability distribution of x is ob¬ 
tained by accumulating the histogram Hn(x) of x over 


all configurations and all chains of size N, given by 

H n {x) = ^2 4,*' (14) 

config. 

here x stands for R e or R g . Note that an angular average 
over all directions has been included in the accumulat¬ 
ing process of the histogram due to spherical symmetry. 
Thus, the normalized histogram of x is given by 

Hn{x) = Hn{x)/ ^2 Hn(x') = A-kCnx 2 Pn{x) ( 15 ) 

x' 

where Cn is the normalization factor such that 
/»00 

Cn / 47 tx 2 Pn{x) = 1 . ( 16 ) 

Jo 

In order to compare the probability distributions of R e 
and R g between various different chain sizes N and bend¬ 
ing energies £&, instead of R e and R g we use the reduced 
end-to-end distance r e = {R 2 /{R 2 )) 1 ' 2 and the reduced 
gyration radius r g = ( R 2 /(R g )) 1 / 2 . Results of the nor¬ 
malized histogram h]\[(r e ) compared to the theoretical 
prediction 


G{r e ) = 4 nr 2 e C N (^^exp (17) 

using {R 2 ) ~ N and Eq. (fl2l) are shown in Fig. [6] Here 
the normalization factor Cn = 1 for all cases. We see the 
nice data collapse for chain sizes IV > 64 and all bending 
energies £b in Fig. [IJi, b and the data are described by 
a universal scaling function G(r e ) for ideal chains. For 
N = 32 the data deviates from the master curve due to 
the finite-size effect. Results of the normalized histogram 
^JV’(j'g) compared to the theoretical prediction 

F(r g ) = 4:Trr 2 C N exp {-b x r~ ad - b 2 r g ) (18) 

using (R 2 ) = k 2 N where k g is a constant and Eq. (fl3l) 
are shown in Fig. [7] Here the normalization factor Cn, 
parameters b\ = a\k~ ad and b 2 = a 2 k g are determined 
numerically depending on the simulation data, and they 
depend slightly on the chain size N. As N —> 00 they will 
tend to asymptotic values (not shown). Near the peak 
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TABLE III: Estimates of root-mean-square end-to-end distance (Re) 1 / 2 , gyration radius (R 2 ) 1 / 2 , bond length 4 = (b 2 ) 1 / 2 , 
persistence lengths, £ p /£b and £ Pt e/£b, using Eqs. 0 and 0, respectively, Flory’s characteristic ratio from Eq. Q, and the 
dimensionless compressibility C g = pksTKr from Eq. (EH for semiflexible chains in a melt for various values of bending energy 
Eb, and for N = 512 and N = 128. 


(a) 


N 

512 

128 

Eb 

0.0 

1.0 

2.0 

3.0 

0.0 

1.0 

2.0 

3.0 

(Rl) 1/2 

72.17 

87.25 

108.95 

132.82 

35.54 

42.96 

53.59 

65.38 

V6 (R 2 ) 1/2 

72.20 

87.12 

108.69 

132.48 

35.51 

42.84 

53.20 

64.41 

4 

2.63 

2.62 

2.62 

2.61 

2.63 

2.62 

2.62 

2.61 

£ p /£b 

_ 

0.96 

1.64 

2.53 

_ 

0.96 

1.64 

2.53 

£ p ,e/£b 

0.45 

0.95 

1.64 

2.53 

0.45 

0.95 

1.64 

2.53 

Coc 

1.24 

2.06 

3.38 

5.12 

1.24 

2.06 

3.38 

5.12 

C 9 

0.25 

0.24 

0.23 

0.23 

0.25 

0.24 

0.23 

0.23 


(b) 


1 

0.8 

2 0.6 

0.2 
0 

0 0.5 1 1.5 2 2.5 




0 0.5 1 1.5 2 2.5 


(R2/<R2 >) 1/2 


r e = (Re/<Re >) 


1/2 


FIG. 6: Normalized probability distributions of r e = ( R 2 /(RI)) 1 / 2 , hjv(r e ), plotted versus r e for polymer chains in a melt. In 
(a) data are for e = 0.0, and for various values of N, as indicated. In (b) data are for N = 512, and for various values of Eb, as 
indicated. The theoretical prediction G(r e ), {Eq. (I17ll l. (solid curve) is also shown for comparison. 


(r g ~ 0.85), we see the finite-size effect is stronger for the 
distribution of r g . The fitting curve F(r g ) with param¬ 
eters b\ = 0.09, b 2 = 2.29, and Cjv = 1.31 determined 
by the least-square fit for N = 512 and £& = 0.0 and 
&i = 0.075, b 2 = 2.20, and C N = 1.13 for N = 512 and 
Eb = 3.0 are shown in Fig.[7ji,b, respectively for compari¬ 
son. The distribution P(R g ) oc Hjg / R) 2 g has its maximum 
value, i.e. P(R g = R g , m ) = max P(R g ), and the corre¬ 
sponding gyration radius R g , m oc R g oc N v . Therefore, 
using Eq. m, the logarithm of the rescaled probability 
distribution as a function of (R g /R g , m ) is given by 


with 


/ 

= A 

A ( R, 


R, 


■g,m 


Rn 


= In 


R g ,m 

Rn 


a 1 = 


a \l h N’' 


P(R g ,m) 

P{ R g) 


ad 


Rn 


R 


9,m 


and ci 2 = 


1 — d 


Ad (i b N h 


R 


g,m 


(19) 


( 20 ) 


where one fitting parameter A is left. Results of 
In {P(R g ^ m )/P(R g )} plotted versus R g /R gym are pre¬ 
sented in Fig. [7t',d. We see that due to the finite-size 
effect, data for fully flexible chains in a melt shown in 
Fig. 0, d only start to converge for N > 128. As N in¬ 
creases, the systematical errors need to be taken into ac¬ 
count. Using the least square fit, it gives A = 0.99(3). It 
is obvious that the distribution can only be well described 
by Eq. (fT9l) for R g > R n rn,, similar as that was found for 
Gaussian chains in Ref. [231 ] . For polymer melts of a fixed 
chain length N = 512, we see that for R g < R g , m the dis¬ 
tribution remains the same for different stiffnesses, while 
for R g > R g ,m the data for £& = 0 are slightly deviated 
from the data for £{, > 0. We plot the same distribution, 
Eq. (fl9l) with A = 0.99 for comparison. 
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0.5 0.75 1 1.25 1.5 1.75 2 

r g = (Rg/<Rg >) V2 


(c) 



(b) 



0.4 0.7 1 1.3 1.6 1.9 2.2 2.5 2.8 3.1 


FIG. 7: (a)(b) Normalized probability distributions of r g = (R 2 /{R 2 )) 1 ^ 2 , h,N(r g ), plotted versus r g for polymer chains in a 
melt. (c)(d) Logarithm of the rescaled probability distribution of gyration radius, In {P{R g ,m)/P(R g )) as a function of R g /R g ,m- 
In (a)(c) data are for £(, = 0.0, and for various values of N, as indicated. The theoretical prediction F(r g ), Eq. 11181) . (solid 
curve) is also shown for comparison. In (b)(d) data are for N = 512, and for various values of Eh , as indicated. The theoretical 
prediction F(r g ), Eq. (1191) . with A = 0.99 (solid curve) is also shown for comparison. 


D. Structure factor and compressibility 

The collective static structure factor 5(g) of polymer 
melts is defined by the total scattering from the center of 
all monomers inside the box, regardless of whether they 
are linked along a polymer chain or not, 

5(<?) = Jj— ex P[*9 • (ri - rj)} \ . (21) 

tot \»=i j=i / 


lar correlations (5i n ter), thus one can write 
5(g) = 5i n tra(g) T 5i n t e r(g) 


with 


N 


^intra 




N 

\n= 1 i,j=l 


and 


(23) 


5 c (g) 

(24) 


Here N to t = Nn c are the total number of monomers and 
(• • •) represents the average over all independent config¬ 
urations, and over all vectors q of the same size q =\ q\. 
Note that since a simple cubic lattice of size L 3 with 
periodic boundary condition is considered for our simu¬ 
lations, only the following q are allowed, 

27r 

g = —(ni,n 2 ,n 3 ) (22) 

where n,; = 0, ±1,±2,... for i = 1, 2, and 3, and 
g > 27 t/L. 5(g) characterizes a competition between the 
intramolecular fluctuations (5i n t ra ) and the intermolecu- 


.. N I n c n c 

5inter = Jf~ ^ ^ ~ *?')] 

iV tOt • • . \ 1 /II 1 

i,J =1 \n=l,n^n' n—1 

(25) 

where the contributions from the intramolecular fluctua¬ 
tions are simply equivalent to the average of the standard 
static structure factor of single polymer chains in a melt, 

5 c (g)- 

The collective structure factors 5(g) which represent 
the conformations of the whole polymer system at dif¬ 
ferent stages from weakly-interactive NRRW chains to 
equilibrated SAW chains in a melt for N = 512 and for 
£b = 0.0, 1.0, 2.0, and 3.0 are shown in Fig. [HI Results 
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(a) (b) 




(c) (d) 




FIG. 8: Structure factor S(q) from the scattering of polymer melts plotted vs. q on log-log scales for eh = 0.0 (a), 1.0 (b), 2.0 
(c) and 3.0 (d). Results are for the conformations of polymer melts at different stages, initial NRRWs (NRRWs), rearranged 
NRRWs after applying pre-packing process (Pre-packing), SAWs obtained from NRRWs by the push-off process before applying 
the pre-packing procedure (EV), and after (Pre-packing+EV), and equilibrated SAWs after polymer melts reaching equilibrium 
(Equilibrium). The structure factor S c (q) for single NRRWs chains generated initially and the Scaling predictions for Gaussian 
coils, S(q) ~ q 2 , and for rigid rods, S(q) ~ q -1 , are also shown for comparison. 


include the initial NRRW chains (NRRWs), rearranged 
NRRW chains through the pre-packing process (Pre¬ 
packing), initial SAW chains obtained from the initial 
NRRW chains by switching on the excluded volume effect 
before (EV) and after (Pre-packing+EV) applying the 
pre-packing procedure, and the equilibrated SAW chains 
in a melt (Equilibrium). Each curve shows the result 
averaged over 32 independent configurations at the in¬ 
termediate states, while 1500 independent configurations 
are considered for the measurement in equilibrium. The 
average static structure factor of single NRRW chains 
generated initially, S c (q), is also included for comparison. 
In Fig. [8] we see that in all cases there is no interactions 
between NRRWs generated in the box at the beginning, 
so the structure factor shows the same scaling behavior 
for single NRRWs in dilute solution, i.e. S(q) = S c (q). 
The scaling predictions S(q) ~ q~ 1 / Vid with = 1/2 for 
ideal chains and S(q) ~ g _1 for rigid rods are verified for 
fully flexible and moderately stiff chains. After the pre¬ 
packing process where NRRW chains are rearranged, the 
collective structure factor decreases in the low-q regime. 
Once all overlapping monomers are pushed off there is no 


difference between the configurations generated by the 
push-off process with and without the pre-packing in the 
intermediate and large q, but for very small q a small 
discrepancy exists between those data sets. It is consis¬ 
tent with our previous observation of ( R 2 (s )) in Fig. [3] 
that the discrepancy becomes more prominent at larger 
scales. 

In the following we only focus on polymer melts con¬ 
sisting of n c = 256 chains of N = 512 monomers in equi¬ 
librium. Figure [9] shows the results for the scattering 
from polymer melts in equilibrium. Four different choices 
of the stiffness characterized by the bending energy £& are 
included. At small q there is a systematic dependence on 
Eb, reflecting the (rather weak) dependence of the com¬ 
pressibility on £{,. It is remarkable that the structure fac¬ 
tor of large values of q is completely independent of the 
chain stiffness. As q increases, the first peak, the so-called 
amorphous halo for non-crystalline materials, appears at 
q ss 7r measuring the mean inter-particle distance (~ 2 
lattice spacings) in a polymer melt. The second peak, 
the so-called Bragg peak, appears at q « 27t probing 
the structure on the scale of the monomers (~ 1 lattice 
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FIG. 9: Structure factors S(q) of polymer melts, plotted vs. 
q on semi-log scales Data are for polymer melts consisting of 
n c = 256 chains of lengths N = 512, and for various values of 
Eb, as indicated. In (a) the peak appears at q ~ 7r measuring 
the mean inter-particle distance in the polymer melts, and 
the second peak aX q ~ 2n is the so-called Bragg peak. 


spacing). The scattering from single chains in a melt in 
equilibrium is shown in Fig. [TOJ In Fig. [TOlone sees that 
S c (q) ~ N ex.p(—q 2 (R 2 )/3) « iV(l — q 2 (R 2 )/ 3) for small 
q in the Guinier regime, then a cross-over occurs to the 
power law of Gaussian coils (ideal chains), S(q) ~ q -1 /" id 
with i/;d = 1/2. For moderately stiff chains one observes 
also a rigid-rod regime S(q) ~ q~ l . The two peaks at 
q ss 7r and q ss 27 t show up in a similar way as that for 
the collective structure S(q) in Fig. [2] and are due to the 
underlying lattice model. In order to clarify whether sin¬ 
gle chains in a melt behave as ideal chains we show the 
structure factors S c (q) in a Kratky-plot in Fig. fToT bl. 
The Debye function [§, 36l-l38| describing the scattering 
from Gaussian chains, 

c c '> o 7 ?- 1 + ex P(“ 7 ?) VI 2 / r>2\ 

C>Debye(<?) = 2-j- With V = Q {Rg} , 

. V . (26) 

is also presented in Fig. 110b for comparison. For fully 

flexible chains (Eb = 0.0) in a melt we see the discrep¬ 
ancy between our simulation results and the Debye func¬ 
tion at the intermediate values of q , which agrees with 
the previous finding [20( | that polymer chains in a melt 
are not random walks. While for Eb = 0.0 the deviation 
from ideality is clearly recognized due to the minimum 
in the Kratky-plot near q = 1, for Eb > 1.0 there oc¬ 
curs no minimum any longer (for N = 512). While the 
deviations from Gaussian statistics are still obvious from 
(cos 0(s)), Fig.[5j they are not easy to extract from S c (q). 
For Eb = 1.0 we see the Kratky plateau over all interme¬ 
diate regimes, while the regime decreases as Eb increases 
since the rigid-rod behavior takes over. 

The isothermal compressibility of kt characterizing 
the degree of density fluctuations on large scale is a mea¬ 
sure for the response of a s yste ms volume to an increase 
in the system pressure [H,j4C|, kt = —(1 /V)(dV/dp)T, 


and is related to the structure factor S(q) as q —> 0, 

lim S(q) = pksTKT = C g (27) 


where C g is defined as the “dimensionless compress¬ 
ibility”. In Fig. [9] we see that as q —> 0 the es¬ 
timates of the collective structure factor for polymer 
chains of size N = 512 in a melt reach a plateau value 
liniq-j-o S(q) ss 0.25. Due to the finite-size effect and lat¬ 
tice artifact (q m i n = 2tt/L), we take the same data but 
plot 1 /S(q) versus q 2 in Fig. [TIJ We see the fluctua¬ 
tions of data near q 2 —> 0. However, the best fit of the 
straight line going through our data gives the estimate of 
C g = [lim^^g l/5(g)] . Results are listed in Table Hill 

We see that C g £ [0.23, 0.25] and depends only weakly 
on the stiffness of the chains. 

In the thermodynamic limit by taking the limit N —> 
00 (but keeping the density p = N/V fixed), then taking 
the limit q —> 0, the local density fluctuations in subvo- 
lumina I^ub < V [H, Hoj 


(n(a) 2 ) - (n(a)) 2 
(n(cr)) 


k B Tpn T = C g 


(28) 


where K su b = 47rcr 3 /3. Choosing cr = 12 such that q m - ln ~ 

_1 I q 

Vsub , i.e. in the limit q —I g m i n = 27 r/L, the local 
density fluctuations estimated from the configurations of 
polymer melts in equilibrium are listed in Table |I| The es¬ 
timates of C g £ [0.27, 0.28] are a little bit higher than the 
estimates from the collective structure factor, but they 
are compatible. We also check the finite-size effect on the 
estimate of the local density fluctuations. Fig. [T2] shows 
the local density fluctuations [(n(cr) 2 ) — (n(a)) 2 ]/(n(a)) 
as a function of a. The arrow and dashed lines shown 
in Fig. [T3] indicate the estimate of the compressibility 
c g (L,Eb) at cr = cr m j n « (47r/3) -1 / 3 (27r/L) -1 . We see 
that the local density fluctuations at certain values of a 
(a < 30) do not depend on the chain lengths N for chains 
of different stiffnesses, while the discrepancy between the 
data sets decreases as the stiffness of chains increases (see 
Fig. fT2hh For cr < cr min the local density fluctuations 
tend to diverge. As the size of polymer melts increases 
(L increases, but p = N/L 3 is fixed) the curves toward 
to a flat curve for cr > cr m j n (see Fig. m>). According to 
the finite-size scaling the estimate of c g can be obtained 
by extrapolating the data of c g (L,£b ) to L _1 —> 0. 


IV. CONCLUSIONS 

We have studied the conformations of polymer chains 
consisting of fully flexible and semiflexible chains in a 
melt based on the bond fluctuation model at volume frac¬ 
tion <f> = 0.5. The initial configurations of polymer melts 
are prepared through several procedures. We first gen¬ 
erate NRRWs in a box with periodic boundary condi¬ 
tions in all three directions, rearrange NRRWs by the 
pre-packing process to reduce the local density fluctua¬ 
tions or remain NRRWs at the same positions, and then 













13 


(a) (b) 




FIG. 10: (a) Structure factors of single chains in a melt, S c (q), plotted vs. q on log-log scales, (b) Same data as in (a) but in 
a Kratky-Plot. Data are for polymer melts consisting of n c = 256 chains of N = 512 monomers, and for various values of Eb, 
as indicated. In (a) the theoretical predictions S c (q) = N(1 — q 2 {Rg)/ 3) at the Guinier regime for small q, S c (q) ~ q~ 2 for a 
Gaussian coil, and S c (q) ~ g _1 for a rigid rod are shown by solid curves for comparison. In (b) the Debye function, Eq. (12611 . 
is also shown for comparison. 


(a) 



(b) 



FIG. 11: Inverse of structure factor of polymer melts, 1/S(q), plotted vs. q 2 for Eb = 0.0 (a), and Eb = 1.0, 3.0 (b). Several 
chain lengths (N — 1) are chosen, as indicated. The asymptotic behavior of 1/S(q) as q 2 m 0 is described by straight lines. 
Values of the dimensionless compressibility C g = [lim 9 2 ^ 0 ys{q)\ 1 are listed in Table Hill 


push off overlapping monomers blocking the same lattice 
sites. We compare the conformations of polymer melts 
at the end of each process to the result obtained from 
the equilibrated configurations. Applying the additional 
pre-packing process seems to have no significant effect 
on preparing initial configurations of polymer chains in a 
melt based on the lattice model after the excluded volume 
interactions are switched on completely. Namely, the es¬ 
timates of the mean square internal end-to-end distance 
(. R 2 (s )) and the collective structure factors S(q) from the 
configurations generated by the two methods (denoted by 
Pre-packing+EV, and EV in Figs. [3] and [5]) are almost 
the same, and very close to estimates from the equili¬ 
brated configurations. There is also no difference of the 
equilibrating time starting from the initial configurations 
generated by these two methods. If, however, one accepts 
a marginally incomplete elimination of excluded volume 


violations (N over « 100 for a system of n c N = 131072 
monomers) one would have significant advantages. An¬ 
other possibility of testing this pre-packing process for 
the future work would be to use it as a criterion of putting 
chains into the box at the beginning. 

In our simulations we combine the algorithms of the 
local 26 moves, slithering-snake, and pivot moves (in¬ 
stead of double-bridging moves |lMl25l) for equilibrating 
the system of polymer melts. Although the total num¬ 
ber of monomers Wot = 131072 which is eight times 
smaller than Wot = 1048576 in Ref. [25| based on the 
same model for £b = 0, the estimates of the mean square 
end-to-end distance (i? 2 ), the mean square gyration ra¬ 
dius (Rg), the mean square bond lengths (b 2 ), and the 
dimensionless compressibility C g are all in perfect agree¬ 
ment with those results given in Ref. [25[ at fixed chain 
size N. For fully flexible and moderately stiff chains the 
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(a) (b) 




FIG. 12: Local density fluctuations [(n(cr) 2 3 ) — (n(cr)} 2 ] /(ra(cr)) plotted as a function of a for polymer chains of different 
stiffnesses in a melt (a) and for fully flexible chains of fixed chain length N = 128 in a melt (b). Three values of chain length N 
and four values of bending energy £(, are chosen, as indicated in (a). Three different sizes of lattice box are chosen, as indicated 
in (b). The locations of a = cr m m = (47r/3) _1,/3 (27r/.L) _1 for various values of L are indicated by an arrow in (a) and dashed 
lines in (b). 


ratio {R 2 )/{R 2 ) « 6 as expected for ideal chains. For 
moderately stiff chains in a melt the internal mean square 
end-to-end distance (R 2 (s)) is well described by the freely 
rotating chain model. Results of the probability distribu¬ 
tions of reduced end-to-end distance r e = (i? 2 /^ 2 )) 1 / 2 
and reduced gyration radius r g = (R 2 g /{R 2 g )) 1 ^ 2 for poly¬ 
mer chains in a melt for various values of N and £& show 
the nice data collapse, and are described by universal 
functions, Eqs. G3 and m , for ideal chains. The col¬ 
lective structure factors S(q) for the whole polymer melts 
and the standard structure factor S c (q) for single chains 
in a melt are also calculated and compared with the theo¬ 
retical predictions. A detailed investigation of S c (q) in a 
Kratky-plot for £5 = 0.0, 1.0, 2.0, and 3.0 shows that for 
fully flexible chains in a melt there is a significant devia¬ 
tion from the Debye function for Gaussian chains at the 
intermediate values of q as found by Wittmer et al. 0 , 
while for £& = 1.0 it is perfectly described by the De¬ 
bye function. Since real polymers (such as polystyrene) 
exhibit some local chain stiffness, it is clear that the de¬ 
viations from Gaussian statistics found for fully flexible 
chains in melts in the Kratky-plot will be very difficult 
to test experimentally. However, since only chains up to 
N = 512 are considered in our simulations, such a devia¬ 
tion may also occur as the chain length increases. Careful 
investigation of the structure factor and the mean square 
internal end-to-end distance for much longer semiflexi- 


ble polymer chains in a melt will be required to clarify 
this matter. The dimensionless compressibility C g de¬ 
termined by the estimates of the collective structure fac¬ 
tors S(q) for small q and the local density fluctuations 
((n 2 ) — (n) 2 )/(n) within a sphere of radius a are com¬ 
patible from our simulations for all cases. We hope that 
the present work will help to the further development 
of a coarse-graining approach by using the BFM as an 
underlying microscopic model and will be useful for the 
interpretation of corresponding experiments searching for 
excluded volume effect in the scattering function of poly¬ 
mers in melts. 
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